Table of Contents

1  Read and process data
2  Align date ranges

In [1]:
import pandas as pd
import numpy as np
import os
import datetime
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from sklearn import tree

import pytz
import itertools
import visualize
import utils
import pydotplus
import xgboost as xgb

from sklearn import metrics

import pvlib
import cs_detection

import visualize_plotly as visualize

from IPython.display import Image

%load_ext autoreload
%autoreload 2

np.set_printoptions(precision=4)
%matplotlib notebook

Read and process data

Only making ground predictions using PVLib clearsky model and statistical model. NSRDB model won't be available to ground measurements.

In [2]:
nsrdb = cs_detection.ClearskyDetection.read_pickle('srrl_nsrdb_1.pkl.gz')
nsrdb.df.index = nsrdb.df.index.tz_convert('MST')
In [3]:
nsrdb.time_from_solar_noon('Clearsky GHI pvlib', 'tfn')
In [4]:
ground = cs_detection.ClearskyDetection.read_pickle('srrl_ground_1.pkl.gz')
ground.df.index = ground.df.index.tz_convert('MST')
ground.time_from_solar_noon('Clearsky GHI pvlib', 'tfn')
In [5]:
feature_cols = [
'tfn',
'abs_ideal_ratio_diff',
'abs_ideal_ratio_diff mean',
'abs_ideal_ratio_diff std',
'abs_ideal_ratio_diff max',
'abs_ideal_ratio_diff min',
'GHI Clearsky GHI pvlib gradient ratio', 
'GHI Clearsky GHI pvlib gradient ratio mean', 
'GHI Clearsky GHI pvlib gradient ratio std', 
'GHI Clearsky GHI pvlib gradient ratio min', 
'GHI Clearsky GHI pvlib gradient ratio max', 
'GHI Clearsky GHI pvlib gradient second ratio', 
'GHI Clearsky GHI pvlib gradient second ratio mean', 
'GHI Clearsky GHI pvlib gradient second ratio std', 
'GHI Clearsky GHI pvlib gradient second ratio min', 
'GHI Clearsky GHI pvlib gradient second ratio max', 
'GHI Clearsky GHI pvlib line length ratio',
'GHI Clearsky GHI pvlib line length ratio gradient',
'GHI Clearsky GHI pvlib line length ratio gradient second'
]

target_cols = ['sky_status']

Align date ranges

In [6]:
ground.df.index[0], ground.df.index[-1]
Out[6]:
(Timestamp('2001-08-04 00:00:00-0700', tz='MST', freq='T'),
 Timestamp('2014-12-31 23:59:00-0700', tz='MST', freq='T'))
In [7]:
nsrdb.df.index[0], nsrdb.df.index[-1]
Out[7]:
(Timestamp('1998-01-01 00:00:00-0700', tz='MST', freq='30T'),
 Timestamp('2015-12-31 23:30:00-0700', tz='MST', freq='30T'))
In [8]:
ground2 = cs_detection.ClearskyDetection(ground.df)
In [9]:
ground2.trim_dates('01-01-2008', '01-01-2012')
ground2.df = ground2.df[ground2.df.index.minute % 30 == 0]
In [10]:
nsrdb2 = cs_detection.ClearskyDetection(nsrdb.df)
In [11]:
nsrdb2.trim_dates('01-01-2008', '01-01-2012')
nsrdb2.df = nsrdb2.df[nsrdb2.df.index.minute % 30 == 0]
In [12]:
vis = visualize.Visualizer()
vis.add_line_ser(ground2.df['GHI'], 'Grnd GHI')
vis.add_line_ser(ground2.df['Clearsky GHI pvlib'], 'Grnd GHIcs')
vis.add_line_ser(nsrdb2.df['GHI'], 'NSRDB GHI')
vis.add_line_ser(nsrdb2.df['Clearsky GHI pvlib'], 'NSRDB GHIcs')
vis.show()
In [13]:
print(np.sqrt(np.mean((nsrdb2.df['GHI'] - ground2.df['GHI'])**2)))

98.7975011251
In [14]:
print(np.sqrt(np.mean((nsrdb2.df['GHI'] - nsrdb2.df['Clearsky GHI pvlib'])**2)))

168.392831054
In [15]:
print(np.sqrt(np.mean((ground2.df['GHI'] - ground2.df['Clearsky GHI pvlib'])**2)))

175.813854654
In [16]:
ground2 = cs_detection.ClearskyDetection(ground.df)
In [17]:
ground2.trim_dates('01-01-2008', '01-01-2012')
In [18]:
nsrdb2 = cs_detection.ClearskyDetection(nsrdb.df)
In [19]:
nsrdb2.trim_dates('01-01-2008', '01-01-2012')
In [20]:
ground2.df['sky_status'] = nsrdb2.df['sky_status']
In [21]:
ground2.df['sky_status'].head()
Out[21]:
2008-01-01 00:00:00-07:00    False
2008-01-01 00:01:00-07:00      NaN
2008-01-01 00:02:00-07:00      NaN
2008-01-01 00:03:00-07:00      NaN
2008-01-01 00:04:00-07:00      NaN
Freq: T, Name: sky_status, dtype: object
In [22]:
ground2.df['sky_status'] = ground2.df['sky_status'].fillna(False)
In [23]:
ground2.scale_model('GHI', 'Clearsky GHI pvlib', 'sky_status')
In [24]:
ground3 = cs_detection.ClearskyDetection(ground2.df)
ground3.trim_dates('07-01-2011', '07-15-2011')
In [25]:
vis = visualize.Visualizer()
vis.add_line_ser(ground3.df['GHI'], 'Grnd GHI')
vis.add_line_ser(ground3.df['Clearsky GHI pvlib'], 'Grnd GHIcs')
vis.add_circle_ser(ground3.df[ground3.df['sky_status']]['GHI'], 'Clear')
vis.show()
In [26]:
utils.calc_all_window_metrics(ground2.df, 11, 'GHI', 'Clearsky GHI pvlib', overwrite=True)
In [27]:
ground2_train = ground2.df[ground2.df.index.minute % 30 == 0]
In [28]:
ground2_train = ground2_train[ground2_train.index < '07-01-2011']
In [29]:
from sklearn import tree
from sklearn import ensemble
import xgboost as xgb
# clf = ensemble.RandomForestClassifier(n_estimators=100)
# clf = xgb.XGBClassifier(max_depth=3, n_estimators=300, learning_rate=.005, reg_lambda=.01)
clf = ensemble.RandomForestClassifier(max_depth=6, n_estimators=128, n_jobs=-1)
In [30]:
from sklearn import model_selection
In [31]:
scores = model_selection.cross_val_score(clf, ground2_train[feature_cols].values, ground2_train[target_cols].values.flatten())
In [32]:
np.mean(scores), np.std(scores)
Out[32]:
(0.92206669276951192, 0.0062942253079680483)
In [33]:
clf.fit(ground2_train[feature_cols].values, ground2_train[target_cols].values.flatten())
Out[33]:
RandomForestClassifier(bootstrap=True, class_weight=None, criterion='gini',
            max_depth=6, max_features='auto', max_leaf_nodes=None,
            min_impurity_split=1e-07, min_samples_leaf=1,
            min_samples_split=2, min_weight_fraction_leaf=0.0,
            n_estimators=128, n_jobs=-1, oob_score=False,
            random_state=None, verbose=0, warm_start=False)
In [34]:
ground2_test = cs_detection.ClearskyDetection(ground2.df)
In [35]:
ground2_test.trim_dates('07-01-2011', '07-15-2011')
In [36]:
pred = ground2_test.iter_predict_daily(feature_cols, 'GHI', 'Clearsky GHI pvlib', clf, 11)
In [37]:
pred2 = pred[pred.index.minute % 30 == 0]
In [38]:
sky_stat = ground2_test.df[ground2_test.df.index.minute % 30 == 0]['sky_status']
In [39]:
cm = metrics.confusion_matrix(sky_stat, pred2)
In [40]:
vis = visualize.Visualizer()
In [41]:
vis.plot_confusion_matrix(cm, labels=['cloudy', 'clear'])
In [42]:
ghi = ground2_test.df['GHI']
ghi_cs = ground2_test.df['Clearsky GHI pvlib']
is_clear = ground2_test.df['sky_status iter']
sky = ground2_test.df['sky_status']
In [43]:
ghi = ghi[ghi.index < '07-15-2011']
ghi_cs = ghi_cs[ghi_cs.index < '07-15-2011']
sky = sky[sky.index < '07-15-2011']
In [44]:
vis = visualize.Visualizer()
vis.add_line_ser(ghi, 'ghi')
vis.add_line_ser(ghi_cs, 'ghi_cs')
vis.add_circle_ser(ghi[is_clear], 'clear')
vis.add_circle_ser(ghi[sky], 'nsrdb')
# vis.add_circle_ser(ghi[ghi.index.isin(sky.index)], 'clear nsrdb')
vis.show()
In [ ]:
In [ ]: